## R Script Output -------------------------------------------------------------
# Appendix Figure B.1: Parallel Trends and DiD Estimates


## Instructions ----------------------------------------------------------------
# Step 1: Adjust MAIN_DIR to where README.txt is located
# Step 2: Run entire script


## setup -----------------------------------------------------------------------
# clean slate
rm(list = ls())
date()

# load packages
pkg <- c("tidyverse",
         "RColorBrewer", 
         "gridExtra",
         "lfe", 
         "stargazer")

lapply(pkg, require, character.only = TRUE)

# set main directory
MAIN_DIR <- "~/Dropbox/Research/ISQ-frei-replication/"


## load data -------------------------------------------------------------------
load(file = paste(MAIN_DIR, "data-merge-zcta.RData", sep = ""))

# subset 
df <- zcta.df 

# subset
df.12.13 <- df %>%
  filter(year == 2012 | year == 2013)

df.11.12 <- df %>%
  filter(year == 2011 | year == 2012)

df.12.14 <- df %>%
  filter(year == 2012 | year == 2014)

df.12.15 <- df %>%
  filter(year == 2012 | year == 2015)

df.12.16 <- df %>%
  filter(year == 2012 | year == 2016)

df.12.17 <- df %>%
  filter(year == 2012 | year == 2017)

df.12.18 <- df %>%
  filter(year == 2012 | year == 2018)

df.12.19 <- df %>%
  filter(year == 2012 | year == 2019)

df.12.20 <- df %>%
  filter(year == 2012 | year == 2020)


## Figure B.1 (left) -----------------------------------------------------------
# summarize quantities of interest
df.high <- df %>%
  group_by(cn_ug_sqmi_bi_2012, year) %>%
  summarize(ln_zhvi_mean = mean(ln_zhvi_mean, na.rm = TRUE)) %>%
  spread(cn_ug_sqmi_bi_2012, ln_zhvi_mean) %>%
  mutate(shift = `CN-UG-high`[year == 2012] - `CN-UG-low`[year == 2012],
         Counterfactual = `CN-UG-low` + shift) %>%
  rename(`Low Exposure` = `CN-UG-low`,
         `High Exposure` = `CN-UG-high`) %>%
  gather(exposure_group, value, -shift, -year) %>%
  mutate(exposure_group = factor(exposure_group,
                                 levels = c("Low Exposure", 
                                            "Counterfactual", 
                                            "High Exposure")))

# set parameters
linetype.high <- c("dotted", "dotdash", "solid")

axis.title.size <- 14

# plot
p.high <- ggplot(df.high, 
                 aes(x = year, 
                     y = value,
                     group = factor(exposure_group))) +
  geom_line(aes(linetype = as.factor(exposure_group))) + 
  scale_linetype_manual(values = linetype.high) +
  geom_vline(xintercept = 2012.5, 
             color = "red") +
  geom_point() +
  scale_x_continuous("Year",
                     breaks = seq(min(df$year), max(df$year), by = 1),
                     labels = seq(min(df$year), max(df$year), by = 1)) +
  scale_y_continuous("Log Typical Home Value (USD)") +
  theme_bw()  +
  theme(plot.title = element_text(size = axis.title.size + 2,
                                  face = "bold",
                                  margin = margin(0, 0, 20, 0),
                                  hjust = 0.5),
        axis.title.y = element_text(size = axis.title.size + 2,
                                    margin = margin(0, 20, 0, 0)),
        axis.title.x = element_text(size = axis.title.size - 1,
                                    margin = margin(20, 0, 0, 0)),
        axis.text = element_text(size = axis.title.size - 2),
        strip.text = element_text(size = axis.title.size - 1),
        strip.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        legend.position = "none") + 
  annotate("segment", 
           x = 2012.6, xend = 2014.5, 
           y = 12.4, yend = 12.4, 
           colour = "black",
           arrow = arrow(length = unit(0.2,"cm"), 
                         ends = "last", type = "closed")) +
  annotate("text", x = 2012.6, y = 12.3, label = "China's\nAnti-Corruption\nCampaign", hjust = 0) +
  annotate("text", x = 2017.2, y = 11.68, label = "Low Exposure", hjust = 0) + 
  annotate("text", x = 2017.2, y = 12.1, label = "Counterfactual", hjust = 0) +
  annotate("text", x = 2017.2, y = 12.31, label = "High Exposure", hjust = 0) 

# save
pdf(file = paste(MAIN_DIR, "Figure-B1-left.pdf", sep = "/"),
    width = 5.5, height = 4.5)
print(p.high)
dev.off()


## Figure B.1 (right) -----------------------------------------------------------
# fit models to data
# 2012 vs. 2013
dd.f.12.13 <- felm(ln_zhvi_mean ~ 
                     anticorrupt +
                     cn_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     cn_g_sqmi_bi_2012_high:anticorrupt +
                     ind_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     share_pop_zcta_female * anticorrupt +
                     
                     share_pop_5_17_zcta * anticorrupt +
                     share_pop_18_24_zcta * anticorrupt + 
                     share_pop_25_34_zcta * anticorrupt +
                     share_pop_35_44_zcta * anticorrupt +
                     share_pop_45_54_zcta * anticorrupt +
                     share_pop_55_64_zcta * anticorrupt +
                     share_pop_65_zcta * anticorrupt + 
                     
                     share_pop_white_zcta * anticorrupt + 
                     share_pop_black_zcta * anticorrupt + 
                     share_pop_hispanic_zcta * anticorrupt + 
                     share_pop_asian_zcta * anticorrupt + 
                     share_pop_cn_zcta * anticorrupt + 
                     
                     share_foreign_born_zcta * anticorrupt + 
                     share_edu_ba_above_zcta * anticorrupt +
                     share_enroll_zcta * anticorrupt + 
                     
                     ln_pop_density_zcta * anticorrupt + 
                     
                     effective_tax_zcta * anticorrupt + 
                     
                     ipw_county * anticorrupt + 
                     
                     employ_rate_zcta * anticorrupt + 
                     median_household_income_zcta * anticorrupt + 
                     vacancy_zcta * anticorrupt 
                   | zcta | 0 | zcta, 
                   keepX = TRUE,
                   data = df.12.13)
summary(dd.f.12.13)


# 2012 vs. 2014
dd.f.12.14 <- felm(ln_zhvi_mean ~ 
                     anticorrupt +
                     cn_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     cn_g_sqmi_bi_2012_high:anticorrupt +
                     ind_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     share_pop_zcta_female * anticorrupt +
                     
                     share_pop_5_17_zcta * anticorrupt +
                     share_pop_18_24_zcta * anticorrupt + 
                     share_pop_25_34_zcta * anticorrupt +
                     share_pop_35_44_zcta * anticorrupt +
                     share_pop_45_54_zcta * anticorrupt +
                     share_pop_55_64_zcta * anticorrupt +
                     share_pop_65_zcta * anticorrupt + 
                     
                     share_pop_white_zcta * anticorrupt + 
                     share_pop_black_zcta * anticorrupt + 
                     share_pop_hispanic_zcta * anticorrupt + 
                     share_pop_asian_zcta * anticorrupt + 
                     share_pop_cn_zcta * anticorrupt + 
                     
                     share_foreign_born_zcta * anticorrupt + 
                     share_edu_ba_above_zcta * anticorrupt +
                     share_enroll_zcta * anticorrupt + 
                     ln_pop_density_zcta * anticorrupt + 
                     
                     effective_tax_zcta * anticorrupt + 
                     
                     ipw_county * anticorrupt + 
                     
                     employ_rate_zcta * anticorrupt + 
                     median_household_income_zcta * anticorrupt + 
                     vacancy_zcta * anticorrupt 
                   | zcta | 0 | zcta, 
                   keepX = TRUE,
                   data = df.12.14)
summary(dd.f.12.14)


# 2012 vs. 2015
dd.f.12.15 <- felm(ln_zhvi_mean ~ 
                     anticorrupt +
                     cn_ug_sqmi_bi_2012_high:anticorrupt +
                    
                     cn_g_sqmi_bi_2012_high:anticorrupt +
                     ind_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     share_pop_zcta_female * anticorrupt +
                     
                     share_pop_5_17_zcta * anticorrupt +
                     share_pop_18_24_zcta * anticorrupt + 
                     share_pop_25_34_zcta * anticorrupt +
                     share_pop_35_44_zcta * anticorrupt +
                     share_pop_45_54_zcta * anticorrupt +
                     share_pop_55_64_zcta * anticorrupt +
                     share_pop_65_zcta * anticorrupt + 
                     
                     share_pop_white_zcta * anticorrupt + 
                     share_pop_black_zcta * anticorrupt + 
                     share_pop_hispanic_zcta * anticorrupt + 
                     share_pop_asian_zcta * anticorrupt + 
                     share_pop_cn_zcta * anticorrupt + 
                     
                     share_foreign_born_zcta * anticorrupt + 
                     share_edu_ba_above_zcta * anticorrupt +
                     share_enroll_zcta * anticorrupt + 
                     ln_pop_density_zcta * anticorrupt + 
                     
                     effective_tax_zcta * anticorrupt + 
                     
                     ipw_county * anticorrupt + 
                     
                     employ_rate_zcta * anticorrupt + 
                     median_household_income_zcta * anticorrupt + 
                     vacancy_zcta * anticorrupt 
                   | zcta | 0 | zcta, 
                   keepX = TRUE,
                   data = df.12.15)
summary(dd.f.12.15)


# 2012 vs. 2016
dd.f.12.16 <- felm(ln_zhvi_mean ~ 
                     anticorrupt +
                     cn_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     cn_g_sqmi_bi_2012_high:anticorrupt +
                     ind_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     share_pop_zcta_female * anticorrupt +
                     
                     share_pop_5_17_zcta * anticorrupt +
                     share_pop_18_24_zcta * anticorrupt + 
                     share_pop_25_34_zcta * anticorrupt +
                     share_pop_35_44_zcta * anticorrupt +
                     share_pop_45_54_zcta * anticorrupt +
                     share_pop_55_64_zcta * anticorrupt +
                     share_pop_65_zcta * anticorrupt + 
                     
                     share_pop_white_zcta * anticorrupt + 
                     share_pop_black_zcta * anticorrupt + 
                     share_pop_hispanic_zcta * anticorrupt + 
                     share_pop_asian_zcta * anticorrupt + 
                     share_pop_cn_zcta * anticorrupt + 
                     
                     share_foreign_born_zcta * anticorrupt + 
                     share_edu_ba_above_zcta * anticorrupt +
                     share_enroll_zcta * anticorrupt + 
                     
                     ln_pop_density_zcta * anticorrupt + 
                     
                     effective_tax_zcta * anticorrupt + 
                     
                     ipw_county * anticorrupt + 
                     
                     employ_rate_zcta * anticorrupt + 
                     median_household_income_zcta * anticorrupt + 
                     vacancy_zcta * anticorrupt 
                   | zcta | 0 | zcta,
                   keepX = TRUE,
                   data = df.12.16)
summary(dd.f.12.16)


# 2012 vs. 2017
dd.f.12.17 <- felm(ln_zhvi_mean ~ 
                     anticorrupt +
                     cn_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     cn_g_sqmi_bi_2012_high:anticorrupt +
                     ind_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     share_pop_zcta_female * anticorrupt +
                     
                     share_pop_5_17_zcta * anticorrupt +
                     share_pop_18_24_zcta * anticorrupt + 
                     share_pop_25_34_zcta * anticorrupt +
                     share_pop_35_44_zcta * anticorrupt +
                     share_pop_45_54_zcta * anticorrupt +
                     share_pop_55_64_zcta * anticorrupt +
                     share_pop_65_zcta * anticorrupt + 
                     
                     share_pop_white_zcta * anticorrupt + 
                     share_pop_black_zcta * anticorrupt + 
                     share_pop_hispanic_zcta * anticorrupt + 
                     share_pop_asian_zcta * anticorrupt + 
                     share_pop_cn_zcta * anticorrupt + 
                     
                     share_foreign_born_zcta * anticorrupt + 
                     share_edu_ba_above_zcta * anticorrupt +
                     share_enroll_zcta * anticorrupt + 
                     ln_pop_density_zcta * anticorrupt + 
                     
                     effective_tax_zcta * anticorrupt + 
                     
                     ipw_county * anticorrupt + 
                     
                     employ_rate_zcta * anticorrupt + 
                     median_household_income_zcta * anticorrupt + 
                     vacancy_zcta * anticorrupt 
                   | zcta | 0 | zcta,
                   keepX = TRUE,
                   data = df.12.17)
summary(dd.f.12.17)


# 2012 vs. 2018
dd.f.12.18 <- felm(ln_zhvi_mean ~ 
                     anticorrupt +
                     cn_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     cn_g_sqmi_bi_2012_high:anticorrupt +
                     ind_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     share_pop_zcta_female * anticorrupt +
                     
                     share_pop_5_17_zcta * anticorrupt +
                     share_pop_18_24_zcta * anticorrupt + 
                     share_pop_25_34_zcta * anticorrupt +
                     share_pop_35_44_zcta * anticorrupt +
                     share_pop_45_54_zcta * anticorrupt +
                     share_pop_55_64_zcta * anticorrupt +
                     share_pop_65_zcta * anticorrupt + 
                     
                     share_pop_white_zcta * anticorrupt + 
                     share_pop_black_zcta * anticorrupt + 
                     share_pop_hispanic_zcta * anticorrupt + 
                     share_pop_asian_zcta * anticorrupt + 
                     share_pop_cn_zcta * anticorrupt + 
                     
                     share_foreign_born_zcta * anticorrupt + 
                     share_edu_ba_above_zcta * anticorrupt +
                     share_enroll_zcta * anticorrupt + 
                     ln_pop_density_zcta * anticorrupt + 
                     
                     effective_tax_zcta * anticorrupt + 
                     
                     ipw_county * anticorrupt + 
                     
                     employ_rate_zcta * anticorrupt + 
                     median_household_income_zcta * anticorrupt + 
                     vacancy_zcta * anticorrupt 
                   | zcta | 0 | zcta,
                   keepX = TRUE,
                   data = df.12.18)
summary(dd.f.12.18)

# 2012 vs. 2019
dd.f.12.19 <- felm(ln_zhvi_mean ~ 
                     anticorrupt +
                     cn_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     cn_g_sqmi_bi_2012_high:anticorrupt +
                     ind_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     share_pop_zcta_female * anticorrupt +
                     
                     share_pop_5_17_zcta * anticorrupt +
                     share_pop_18_24_zcta * anticorrupt + 
                     share_pop_25_34_zcta * anticorrupt +
                     share_pop_35_44_zcta * anticorrupt +
                     share_pop_45_54_zcta * anticorrupt +
                     share_pop_55_64_zcta * anticorrupt +
                     share_pop_65_zcta * anticorrupt + 
                     
                     share_pop_white_zcta * anticorrupt + 
                     share_pop_black_zcta * anticorrupt + 
                     share_pop_hispanic_zcta * anticorrupt + 
                     share_pop_asian_zcta * anticorrupt + 
                     share_pop_cn_zcta * anticorrupt + 
                     
                     share_foreign_born_zcta * anticorrupt + 
                     share_edu_ba_above_zcta * anticorrupt +
                     share_enroll_zcta * anticorrupt + 
                     ln_pop_density_zcta * anticorrupt + 
                     
                     effective_tax_zcta * anticorrupt + 
                     
                     ipw_county * anticorrupt + 
                     
                     employ_rate_zcta * anticorrupt + 
                     median_household_income_zcta * anticorrupt + 
                     vacancy_zcta * anticorrupt 
                   | zcta | 0 | zcta,
                   keepX = TRUE,
                   data = df.12.19)
summary(dd.f.12.19)

# 2012 vs. 2020
dd.f.12.20 <- felm(ln_zhvi_mean ~ 
                     anticorrupt +
                     cn_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     cn_g_sqmi_bi_2012_high:anticorrupt +
                     ind_ug_sqmi_bi_2012_high:anticorrupt +
                     
                     share_pop_zcta_female * anticorrupt +
                     
                     share_pop_5_17_zcta * anticorrupt +
                     share_pop_18_24_zcta * anticorrupt + 
                     share_pop_25_34_zcta * anticorrupt +
                     share_pop_35_44_zcta * anticorrupt +
                     share_pop_45_54_zcta * anticorrupt +
                     share_pop_55_64_zcta * anticorrupt +
                     share_pop_65_zcta * anticorrupt + 
                     
                     share_pop_white_zcta * anticorrupt + 
                     share_pop_black_zcta * anticorrupt + 
                     share_pop_hispanic_zcta * anticorrupt + 
                     share_pop_asian_zcta * anticorrupt + 
                     share_pop_cn_zcta * anticorrupt + 
                     
                     share_foreign_born_zcta * anticorrupt + 
                     share_edu_ba_above_zcta * anticorrupt +
                     share_enroll_zcta * anticorrupt + 
                     ln_pop_density_zcta * anticorrupt + 
                     
                     effective_tax_zcta * anticorrupt + 
                     
                     ipw_county * anticorrupt + 
                     
                     employ_rate_zcta * anticorrupt + 
                     median_household_income_zcta * anticorrupt + 
                     vacancy_zcta * anticorrupt 
                   | zcta | 0 | zcta,
                   keepX = TRUE,
                   data = df.12.20)
summary(dd.f.12.20)


# 2011 vs. 2012
dd.f.11.12 <- felm(ln_zhvi_mean ~
                     placebo_anticorrupt_2012 +
                     cn_ug_sqmi_bi_2012_high:placebo_anticorrupt_2012 +
                     
                     cn_g_sqmi_bi_2012_high:placebo_anticorrupt_2012 +
                     ind_ug_sqmi_bi_2012_high:placebo_anticorrupt_2012 +
                     
                     share_pop_zcta_female * placebo_anticorrupt_2012 +
                     
                     share_pop_5_17_zcta * placebo_anticorrupt_2012 +
                     share_pop_18_24_zcta * placebo_anticorrupt_2012 + 
                     share_pop_25_34_zcta * placebo_anticorrupt_2012 +
                     share_pop_35_44_zcta * placebo_anticorrupt_2012 +
                     share_pop_45_54_zcta * placebo_anticorrupt_2012 +
                     share_pop_55_64_zcta * placebo_anticorrupt_2012 +
                     share_pop_65_zcta * placebo_anticorrupt_2012 + 
                     
                     share_pop_white_zcta * placebo_anticorrupt_2012 + 
                     share_pop_black_zcta * placebo_anticorrupt_2012 + 
                     share_pop_hispanic_zcta * placebo_anticorrupt_2012 + 
                     share_pop_asian_zcta * placebo_anticorrupt_2012 + 
                     share_pop_cn_zcta * placebo_anticorrupt_2012 + 
                     
                     share_foreign_born_zcta * placebo_anticorrupt_2012 + 
                     share_edu_ba_above_zcta * placebo_anticorrupt_2012 +
                     share_enroll_zcta * placebo_anticorrupt_2012 + 
                     
                     ln_pop_density_zcta * placebo_anticorrupt_2012 + 
                     
                     effective_tax_zcta * placebo_anticorrupt_2012 + 
                     
                     ipw_county * placebo_anticorrupt_2012 + 
                     
                     employ_rate_zcta * placebo_anticorrupt_2012 + 
                     median_household_income_zcta * placebo_anticorrupt_2012 + 
                     vacancy_zcta * placebo_anticorrupt_2012 
                   | zcta | 0 | zcta, 
                   keepX = TRUE,
                   data = df.11.12)
summary(dd.f.11.12)

# set var order
var.time <- c("2011-2012", 
              "2012-2013",
              "2012-2014", 
              "2012-2015",
              "2012-2016",
              "2012-2017",
              "2012-2018",
              "2012-2019",
              "2012-2020")

# extract coefs and clustered SEs
coef.df.time <- tibble(var = factor(var.time, 
                                    levels = var.time),
                       coef = c(coef(dd.f.11.12)["placebo_anticorrupt_2012:cn_ug_sqmi_bi_2012_high"],
                                coef(dd.f.12.13)["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                                coef(dd.f.12.14)["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                                coef(dd.f.12.15)["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                                coef(dd.f.12.16)["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                                coef(dd.f.12.17)["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                                coef(dd.f.12.18)["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                                coef(dd.f.12.19)["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                                coef(dd.f.12.20)["anticorrupt:cn_ug_sqmi_bi_2012_high"]),
                       cse = c(dd.f.11.12$cse["placebo_anticorrupt_2012:cn_ug_sqmi_bi_2012_high"],
                               dd.f.12.13$cse["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                               dd.f.12.14$cse["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                               dd.f.12.15$cse["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                               dd.f.12.16$cse["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                               dd.f.12.17$cse["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                               dd.f.12.18$cse["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                               dd.f.12.19$cse["anticorrupt:cn_ug_sqmi_bi_2012_high"],
                               dd.f.12.20$cse["anticorrupt:cn_ug_sqmi_bi_2012_high"]))

# set alpha level
alpha <- 0.05

# compute C.I.
coef.df.time <- coef.df.time %>%
  mutate(multiplier = qnorm(1 - alpha/2),
         lwr = coef - multiplier * cse,
         upr = coef + multiplier * cse)

# add symbol indicator
coef.time.shape <- c(rep("placebo", 1),
                     rep("effect", 8))

coef.df.time <- cbind(coef.df.time, coef.time.shape)

# df for annotations
txt.shift.time <- 0.004
txt.pos.time <- tibble(var = factor(var.time, levels = var.time),
                       coef = c(coef.df.time[1,]$lwr - txt.shift.time / 2,
                                coef.df.time[2,]$upr + txt.shift.time / 2,
                                coef.df.time[3,]$upr + txt.shift.time / 2,
                                coef.df.time[4,]$upr + txt.shift.time / 2,
                                coef.df.time[5,]$upr + txt.shift.time / 2,
                                coef.df.time[6,]$upr + txt.shift.time / 2,
                                coef.df.time[7,]$upr + txt.shift.time / 2,
                                coef.df.time[8,]$upr + txt.shift.time / 2,
                                coef.df.time[9,]$upr + txt.shift.time / 2),
                       lwr = coef.df.time$lwr,
                       upr = coef.df.time$upr)

# set coefficient symbol
shapes <- c(16, 17)

# set size
axis.title.size <- 16

# plot
p.coef.time <- ggplot(data = coef.df.time,
                      aes(x = var, y = coef,
                          ymin = lwr, 
                          ymax = upr,
                          shape = factor(coef.time.shape))) +
  geom_pointrange(size = 0.75) +
  geom_hline(yintercept = 0,
             color = I(hsv(0/12, 7/12, 7/12)), 
             linetype = "dashed") +
  scale_y_continuous("Estimated DiD Effect on Home Value (USD, log)",
                     limits = c(-0.02, 0.041)) +
  scale_shape_manual(values = shapes) + 
  geom_text(data = txt.pos.time, 
            aes(label = var),
            size = axis.title.size - 12.5) +
  theme_bw() +  
  theme(legend.position = "none",
        axis.title.y = element_text(size = axis.title.size,
                                    margin = margin(0, 20, 0, 0)),
        axis.text.y = element_text(size = axis.title.size),
        strip.background = element_blank(),
        panel.border = element_rect(colour = "black"),
        strip.text = element_text(size = axis.title.size + 3),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        plot.title = element_text(size = axis.title.size - 1,
                                  face = "bold",
                                  margin = margin(0, 0, 10, 0),
                                  hjust = 0.5),
        panel.spacing = unit(1.5, "lines"),
        panel.grid = element_blank()) 

# save
pdf(file = paste(MAIN_DIR, "Figure-B1-right.pdf", sep = "/"), 
    width = 8, height = 5)
print(p.coef.time)
dev.off()

